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Abstract 

We show that the computational effort for the numerical solution of fermionic 
quantum systems, occurring e.g., in quantum chemistry, solid state physics, 
field theory in principle grows with less than the square of the particle number 
for problems stated in one space dimension and with less than the cube of 
the particle number for problems stated in three space dimensions. This is 
proven by representation of effective algorithms for fermion systems in the 
framework of the Feynman Path Integral. 
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I. INTRODUCTION 

When Feynman stated the thermodynamical Path Integral for fermions (see [Q] and 
references therein), he was in doubt whether this formula could ever be used in numerical 
applications. However, in the last few years considerable progress has been made in the 
development of Path Integral Monte Carlo algorithms. Especially D. Ceperly et. al. ^ 
proposed many very useful refinements. Nevertheless, the problem of applying Path Integral 
methods to fermion systems remains hard to be solved. The reasons are the so-called fermion 
sign problem and the computational effort, which grows normally proportional to the fourth 
power of the number of particles N. In this letter we present a formally exact method with 
a computational effort growing only with the cube of the particle number. 
We further show, that there are some strong indications, that the computational effort 
should in general grow only with the square of the particle number. For problems stated in 
one space dimension, we give an explicit proof for that. 



II. EFFECTIVE PATH INTEGRAL METHOD 

A quantum body system in thermal equilibrium can be described completely by the 
quantum partition function 

Z = Trexp(-/3(Ho + Hi)) , (1) 

where Hq and Hi are the kinetic and potential energy operator. The Path Integral formalism 
leads to a computational convenient equation to calculate Z. Using Trotter's formula P| Z 
can be approximated by 

Z = Tt (exp(-^Ho) exp(-^Ho) + 0(/?VM^). (2) 

For a system of polarized fermions in d space dimensions this leads to the dNM dimensional 
integral: 
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and the periodicity condition Xi{M + 1) = afj(l). The particle mass is denoted by m. 
A number of useful modifications to (3), which improve the convergence vs. the number of 
timesteps M (see e.g. [^0), are known. Because these modifications do not affect any of 
the following we will not stress them here. 

Usually the Metropolis algorithm |^ is adopted to evaluate The improvements pre- 
sented in this paper are based on a careful analysis of the numerical algorithms. 



Within the Metropolis Monte Carlo procedure a Markov process has to be generated, 
under which the probability distribution is stationary. In our case this condition is fulfilled 
by the following widely used procedure, in which the sampling of the probability distribution 
is done by generating two different types of random moves of the particle coordinates. 
In every microscopic step, all time slices of all particle coordinates are moved separately 
through 

Xi{a) ^ Xi{a) + Ax , (5) 

where Ax is a randomly chosen vector. In every macroscopic step all time slices of every 
particle coordinate are moved at once by constant vector. 

Xi{a) ^ Xi{a) + Ay, a = 1...M. (6) 

In both cases a move is accepted, if the absolute value of the ratio of weight functions, which 
are simply the integrands in (^, is greater than a homogeneous random number. Here it is 
important that the absolute values have to be taken, because the determinant occurring in 
the weight function W{p) becomes sometimes negative. 
Any observable X can then be evaluated through 
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E^=i sign{W{p)) 



where the summation in (|^) runs over all Metropolis steps, X{p) being the value of X in the 
p-th step. 



Albeit the fact, that this method has been applied in some cases (see e.g. |]^), its useful- 
ness is limited by its high computational costs. The main effort in the numerical computation 
of (|^) is the calculation of the determinant. In every complete microscopic motion the de- 
terminant has to be calculated N x M times. The complete algorithm is therefore of order 
< N'^, if standard matrix factorizations , which are of order A^'^, are used. The lower sign 
stands, because we found that the number of iterations necessary to achieve a given preci- 
sion in the numerical solution decreases with an increasing number of particles. This results 
simply from the fact that a larger number of identical particles yields better statistics. 
Now observe, that every change of a microscopic coordinate affects only the i'th row of 
matrix A{a + 1, a) and the i'th column of matrix A{a, a — 1). This simple fact can be used 
to reduce the numerical effort by a factor n. The reason for that is that the updating of an 
LU factorization of a matrix after a row or column exchange can be done by a numerical 
effort proportional to A^^. Descriptions of such algorithms can be found in some textbooks 
on matrix computations We thus claim that the total algorithm is only of order A^^. 
Encouraged by the above result and through the fact that the matrix has some hidden 
symmetries having only 2dN independent parameters instead of A^^, we found for the spe- 
cial case of systems describable in one space dimension {d = 1) a further reduction of the 
numerical effort. 

Absorbing the trivial constant factors in (^) the problem is reduced to calculate the deter- 
minant of matrices of the form 

Bij = exp{-^{xi-yjf). (8) 
Using the linear properties of det we have 
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det{B) = det(M) J] exp(--(x| + y])) . 
We now show, that the determinant of matrices of the special form 

Mi J = exp^XiUj) 



(9) 
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can be calculated with an effort of only A^^ operations. This can most easily be seen using 
the decomposition M = RST of (|TU]). 
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The truncation of these infinite matrices to simple (A^ x A^)-matrices corresponds to a 
truncation of the exponential series at order A^ — 1. For sufficiently large there is obviously 
no problem in doing so and we have 

det(M) = det(i?) det(5) det(T). (12) 

The calculation of det(S') is just trivial and yields a constant for given A^. R and T are 
Vandermonde matrices, for which the calculation of the determinants is extremely simple 
0]. For example for R we have 

det{R)= n (Xi-Xj). (13) 

l<j<i<N 

The complete calculation of det M thus involves A^(A^ — 1) subtractions and A^(A^ — 1) 
multiplications altogether. Again using the fact, that only one coordinate is changed in one 
microscopic step, the algorithm to compute the determinant is of order A^ and the complete 
Monte Carlo algorithm thus of order A^^. 
Thus we are able to rewrite (|]) as 
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III. CONCLUSION 



Our almost trivial looking results based on basic linear algebra contribute a significant 
improvement to the Path Integral method. The decomposition of (|10|) found together with 
the truncation of the indefinite matrices is somewhat unsatisfactory, because the method 
applies in principle only to large particle numbers. Nonetheless our results should be en- 
couraging to continue to look for efficient algorithms to compute the determinant of (PUf). 
Because of the equivalence of the mathematical problem, quite generally our result implies 
that any solution method for quantum systems involving fermions should be bound by a 
computational effort of order A^^. 

Indeed, a completely different method of computing the Path Integral for fermion and 



boson systems was found, which has this property |10,11 
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